      subroutine cmfdp1

      use param

#include <global.h>
      include 'ntbytes.h'
      include 'geom.h'
      include 'xsec.h'
      include 'cmfdp1.h'
      include 'sanm.h'
      include 'sanmjin.h'
      include 'setls.h'
      include 'eigv.h'
      include 'nem.h'

      real(NBF) btm1,bt1,btp1,btm2,bt2,ppm,pp,ppp,sffdm1,sffdm2,del,ss,rhs &
     ,         bcoefl,bcoefr,sffdml,sffdmr,jnet_l,jnet_r

! flux/current correction factor
      do im=1,ng
        do j=1,tnode
          jm1=j-1
          jp1=j+1
          dtill(j,im)=2*beta(1,jm1,im)*beta(1,j,im)/(beta(1,jm1,im)+beta(1,j,im))
          dtilr(j,im)=2*beta(1,j,im)*beta(1,jp1,im)/(beta(1,j,im)+beta(1,jp1,im))
        enddo
        dtill(1,im)=bcm_p1(1)
        dtilr(tnode,im)=bcm_p1(2)
      enddo
      
      do im=1,ng
        do j=1,tnode
          jm1=j-1
          jp1=j+1
          btm1=beta(1,jm1,im)
          bt1=beta(1,j,im)
          btp1=beta(1,jp1,im)
          ppm=psphi(1,jm1,im)
          pp=psphi(1,j,im)
          ppp=psphi(1,jp1,im)

          dhatl(j,im)=-(jnet(1,j,im)+dtill(j,im)*(pp-ppm))/(pp+ppm)   ! current correction factor
          dhatr(j,im)=-(jnet(1,jp1,im)+dtilr(j,im)*(ppp-pp))/(ppp+pp)
#define NEMP3
#ifdef NEMP3
! NEM P3
          sffdm1=(bt1*pp+btm1*ppm)/(bt1+btm1)
          flcor(j,im)=(sflux(1,j,im)-sffdm1)/(pp+ppm)
          if (j.eq.tnode) then
            sffdm1=(btp1*ppp+bt1*pp)/(btp1+bt1)
            flcor(jp1,im)=(sflux(1,jp1,im)-sffdm1)/(pp+ppp)
          endif
#else
          sffdm1=(bt1*avgflx(1,j,im)+btm1*avgflx(1,jm1,im))/(bt1+btm1)  ! surface flux for the FDM
          flcor(j,im)=(sflux(1,j,im)-sffdm1)/(avgflx(1,j,im)+avgflx(1,jm1,im))  ! flux correction factor
#endif
!          btm2=beta(2,jm1,im)
!          bt2=beta(2,j,im)
!          sffdm2=(bt2*avgflx(2,j,im)+btm2*avgflx(2,jm1,im))/(bt2+btm2)
!          flcor(2,j,im)=(sflux(2,j,im)-sffdm2)/(avgflx(2,j,im)+avgflx(2,jm1,im))
        enddo
        if (bcb.eq.1) then
          dhatl(1,im)=-jnet(1,1,im)/sflux(1,1,im)-dtill(1,im)
        endif
        if (bcu.eq.1) then
          dhatr(tnode,im)=-jnet(1,tnode+1,im)/sflux(1,tnode+1,im)+dtilr(tnode,im)
        endif
      enddo
      
      do im=1,ng
        do j=1,tnode
          czf(1,j,im)=-dtill(j,im)+dhatl(j,im)
          czf(2,j,im)=-dtilr(j,im)-dhatr(j,im)
          digf(j,im)=(dtilr(j,im)-dhatr(j,im)+dtill(j,im)+dhatl(j,im))+xsr(j,im)*hz(j)
        enddo
      enddo

      do im=1,ng
         j=1
         del=digf(j,im)
         delinvf(j,im)=1/del
         jb=j
         do j=2,tnode
           deliauf(jb,im)=delinvf(jb,im)*czf(2,jb,im)
           del=digf(j,im)-czf(1,j,im)*deliauf(jb,im)
           delinvf(j,im)=1/del
           jb=j
         enddo
      enddo

! source term
      do im=1,ng
        do j=1,tnode
          ss=0
          do m=1,ng         
            ss=ss+avgflx(1,j,m)*xssm(j,im)%from(m)
          enddo
!          src(1,j,im)=reigv*xchi(j,im)*psif(j)+ss*hz(j)
          src(1,j,im)=reigv*xchi(j,im)*psif(j)+(ss+2*xsr(j,im)*avgflx(2,j,im))*hz(j)
        enddo
!      enddo
!      do im=1,ng
! forward substitution
        j=1   
        rhs=src(1,j,im)
        y(j)=delinvf(j,im)*rhs
        jb=j
        do j=2,tnode
          rhs=src(1,j,im)-czf(1,j,im)*y(jb)
          y(j)=delinvf(j,im)*rhs
          jb=j
        enddo            
! backward substitution
        j=tnode   
        psphi(1,j,im)=y(j)
        jt=j 
        do j=tnode-1,1,-1
          psphi(1,j,im)=y(j)-deliauf(j,im)*psphi(1,jt,im)
          jt=j
        enddo
        do j=1,tnode
          avgflx(1,j,im)=psphi(1,j,im)-2*avgflx(2,j,im)
        enddo                  
      enddo ! for im          

      call dumflx

! Determine the incoming partial current J_in
      do im=1,ng
        do j=1,tnode
          jm1=j-1
          jp1=j+1
          bt1=beta(1,j,im)
          btm1=beta(1,jm1,im)
          btp1=beta(1,jp1,im)
          pp=psphi(1,j,im)
          ppm=psphi(1,jm1,im)
          ppp=psphi(1,jp1,im)
#define NEMP3
#ifdef NEMP3
! NEM P3 (using a pseudo flux)
          sffdml=(btm1*ppm+bt1*pp)/(btm1+bt1)
          sffdmr=(btp1*ppp+bt1*pp)/(btp1+bt1)
          sfluxl(1)=sffdml+flcor(j,im)*(ppm+pp)
          sfluxr(1)=sffdmr+flcor(jp1,im)*(ppp+pp) 
          sfluxl(2)=sflux(2,j,im)
          sfluxr(2)=sflux(2,jp1,im)
          jnet_l=-dtill(j,im)*(pp-ppm)-dhatl(j,im)*(pp+ppm)
          jnet_r=-dtilr(j,im)*(ppp-pp)-dhatr(j,im)*(ppp+pp)

          ajil(1,j,im)=0.25*sfluxl(1)+0.5*jnet_l-0.1875*sfluxl(2)
          ajir(1,j,im)=0.25*sfluxr(1)-0.5*jnet_r-0.1875*sfluxr(2)
          ajil(2,j,im)=0.4375*sfluxl(2)+0.5*jnet(2,j,im)-0.0625*sfluxl(1)
          ajir(2,j,im)=0.4375*sfluxr(2)-0.5*jnet(2,jp1,im)-0.0625*sfluxr(1)       
#else
          sffdml=(btm1*avgflx(1,jm1,im)+bt1*avgflx(1,j,im))/(btm1+bt1)
          sffdmr=(btp1*avgflx(1,jp1,im)+bt1*avgflx(1,j,im))/(btp1+bt1)
          sfluxl(1)=sffdml+flcor(j,im)*(avgflx(1,jm1,im)+avgflx(1,j,im))
          sfluxr(1)=sffdmr+flcor(jp1,im)*(avgflx(1,jp1,im)+avgflx(1,j,im))  
          sfluxl(2)=sflux(2,j,im)
          sfluxr(2)=sflux(2,jp1,im)
          jnet_l=-dtill(j,im)*(pp-ppm)-dhatl(j,im)*(pp+ppm)
          jnet_r=-dtilr(j,im)*(ppp-pp)-dhatr(j,im)*(ppp+pp)

          jin(1,j,im)=0.25*sfluxl(1)+0.5*jnet_l+0.3125*sfluxl(2)
          jin(2,j,im)=0.25*sfluxr(1)-0.5*jnet_r+0.3125*sfluxr(2)
#endif
        enddo
        if (bcb.eq.1) jin(1,1,im)=0; ajil(1,1,im)=0
        if (bcu.eq.1) jin(2,tnode,im)=0; ajir(1,tnode,im)=0  
      enddo
 
      end subroutine
